#===============================================================================
# The following analyses were carried out using R version 4.0.3
#===============================================================================
# Install packages
# install.packages("tidyverse")
# install.packages("lme4")
# install.packages("stargazer")
#-------------------------------------------------------------------------------
#Load packages
#-------------------------------------------------------------------------------
rm(list = ls())
library(tidyverse)
library(lme4)
library(stargazer)
#===============================================================================
# set working directory
setwd(" ")

# read data
mydata<-read.csv("FixedAssets&Corruption.csv")

#===============================================================================
# Table 1: Baseline Model
#===============================================================================
base_1 <- lm(letcs ~ 
               lintensity_investment_employment +
               marketsize + growthrate + lscale + soe_influence +
               mgovthelp + mtaxrate + 
               soe + collective + private + foreign  + 
               lage + lemp + sales_otherprov + govtsales + 
               soesales + relationship + licenses + interaction +
               + factor(province),data=mydata)
#------------
base_2 <- lm(letcs ~ 
               lintensity_investment_employment +
               marketsize + growthrate + lscale + soe_influence +
               mgovthelp + mtaxrate + 
               soe + collective + private + foreign +
               lage + lemp + sales_otherprov + govtsales + 
               soesales + relationship + licenses + interaction + lceopay +
               factor(province),data=mydata)
#------------
base_3 <- lm(letcs ~ 
               lintensity_investment_employment +
               marketsize + growthrate + lscale + soe_influence +
               mgovthelp + mtaxrate + 
               soe + collective + private + foreign +
               lage + lemp + sales_otherprov + govtsales + 
               soesales + relationship + licenses + interaction + gm_govt +
               factor(province),data=mydata)

#-------------------------------------------------------------------------------
# Calculate Clustered Standard Errors
#-------------------------------------------------------------------------------
#install.packages("lmtest")
#install.packages("sandwich")
library(lmtest)
library(sandwich)

CL.vcov<-function(model, cluster){
  #calculate degree of freedom adjustment
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- length(coef(model)) - M + 1
  dfc <- (M/(M-1))*((N-1)/(N-K))
  #calculate the uj's
  uj  <- apply(estfun(model),2, function(x) tapply(x, cluster, sum))
  #use sandwich to get the var-covar matrix
  vcovCL <- dfc*sandwich(model, meat=crossprod(uj)/N)
  return(vcovCL)
}

#-------------------------------------------------------------------------------

base_1_vcovCL <- CL.vcov(base_1, base_1$model[,"factor(province)"])
base_1_CL <- coeftest(base_1, base_1_vcovCL, df=29)

base_2_vcovCL <- CL.vcov(base_2, base_2$model[,"factor(province)"])
base_2_CL <- coeftest(base_2, base_2_vcovCL, df=29)

base_3_vcovCL <- CL.vcov(base_3, base_3$model[,"factor(province)"])
base_3_CL <- coeftest(base_3, base_3_vcovCL, df=29)

stargazer(base_1_CL,base_2_CL,base_3_CL, omit= c("province"),
          type = "html", out = "Table1.doc",
          covariate.labels =c("Fixed-Asset Intensity (log)","Market Size","Growth Rate",
                              "Scale Economies (log)","SOE Dominance", "Government Help","Tax Burden",
                              "SOE","Collective","Private","Foreign","Age (log)",
                              "Employees (log)","Out-of-Province Sales","Proportion of Sales to Government",
                              "Proportion of Sales to SOEs","Years of Relationship",
                              "Licenses","Interaction with Government","CEO Pay (log)",
                              "Government Appointed General Manager"),no.space = TRUE)

#===============================================================================
# Table 2: Mediation Analysis
#===============================================================================

mediate.data <- mydata %>% dplyr::select(province,prov_ind,
                                      w_con_sal4,hhi,
                                      marketsize,growthrate,lscale, soe_influence,
                                      mgovthelp,mtaxrate,
                                      lintensity_investment_employment) %>% 
                                      distinct(prov_ind, .keep_all = T) %>% 
                                      filter(complete.cases(.)==T)

#-------------------------------------------------------------------------------
# CR4
#-------------------------------------------------------------------------------
model.cr4.0 <- lm(w_con_sal4 ~  
                    lintensity_investment_employment + 
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate + factor(province),
                  data = mediate.data)

cr4.0_vcovCL <- CL.vcov(model.cr4.0, model.cr4.0$model[,"factor(province)"])
cr4.0_CL <- coeftest(model.cr4.0, cr4.0_vcovCL,df=29)

#------------
model.cr4.1 <- lm(letcs ~ 
                    lintensity_investment_employment +
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate + 
                    soe + collective + private + foreign +
                    lage + lemp + sales_otherprov + govtsales + 
                    soesales + relationship + licenses + interaction +
                    + factor(province),data=mydata)

cr4.1_vcovCL <- CL.vcov(model.cr4.1, model.cr4.1$model[,"factor(province)"])
cr4.1_CL <- coeftest(model.cr4.1, cr4.1_vcovCL,df=29)

#------------
model.cr4.2 <- lm(letcs ~ 
                    lintensity_investment_employment +
                    w_con_sal4 + 
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate + 
                    soe + collective + private + foreign +
                    lage + lemp + sales_otherprov + govtsales + 
                    soesales + relationship + licenses + interaction +
                    + factor(province),data=mydata)

cr4.2_vcovCL <- CL.vcov(model.cr4.2, model.cr4.2$model[,"factor(province)"])
cr4.2_CL <- coeftest(model.cr4.2, cr4.2_vcovCL, df=29)


#-------------------------------------------------------------------------------
# HHI
#-------------------------------------------------------------------------------
model.hhi.0 <- lm(hhi ~  
                    lintensity_investment_employment + 
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate + factor(province),
                    data = mediate.data)

hhi.0_vcovCL <- CL.vcov(model.hhi.0, model.hhi.0$model[,"factor(province)"])
hhi.0_CL <- coeftest(model.hhi.0, hhi.0_vcovCL,df=29)

#------------
model.hhi.1 <- lm(letcs ~ 
                    lintensity_investment_employment +
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate + 
                    soe + collective + private + foreign +
                    lage + lemp + sales_otherprov + govtsales + 
                    soesales + relationship + licenses + interaction +
                    + factor(province),data=mydata)

hhi.1_vcovCL <- CL.vcov(model.hhi.1, model.hhi.1$model[,"factor(province)"])
hhi.1_CL <- coeftest(model.hhi.1, hhi.1_vcovCL,df=29)

#------------
model.hhi.2 <- lm(letcs ~ 
                    lintensity_investment_employment +
                    hhi + 
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate + 
                    soe + collective + private + foreign +
                    lage + lemp + sales_otherprov + govtsales + 
                    soesales + relationship + licenses + interaction +
                    + factor(province),data=mydata)

hhi.2_vcovCL <- CL.vcov(model.hhi.2, model.hhi.2$model[,"factor(province)"])
hhi.2_CL <- coeftest(model.hhi.2, hhi.2_vcovCL,df=29)


stargazer(cr4.0_CL,cr4.1_CL,cr4.2_CL,hhi.0_CL,hhi.1_CL,hhi.2_CL, omit= c("province"),
          type = "html", out = "Table2.doc",
          covariate.labels =c("Fixed-Asset Intensity (log)", "Market Concentration (CR4)", 
                              "Market Concentration (HHI)","Market Size","Growth Rate",
                              "Scale Economies (log)","SOE Dominance", "Government Help","Tax Burden",
                              "SOE","Collective","Private","Foreign", "Age (log)",
                              "Employees (log)","Out-of-Province Sales","Proportion of Sales to Government",
                              "Proportion of Sales to SOEs","Years of Relationship",
                              "Licenses","Interaction with Government"),no.space = TRUE)


#===============================================================================
# Table 3: Multilevel Model with Interaction Term
#===============================================================================
MLM_LSC <- lmer(letcs ~ 
                  lintensity_investment_employment +
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate  +
                    LSC + gdpper_1999_2003 + pop +
                    soe + collective + private + foreign +
                    lage + lemp + sales_otherprov + govtsales + 
                    soesales + relationship + licenses + interaction+
                    LSC * lintensity_investment_employment + (1|province),
                    data = mydata)

#------------
MLM_court <- lmer(letcs ~ 
                    lintensity_investment_employment +
                      marketsize + growthrate + lscale + soe_influence +
                      mgovthelp + mtaxrate  +
                      TrustCourts + gdpper_1999_2003 + pop +
                      soe + collective + private + foreign +
                      lage + lemp + sales_otherprov + govtsales + 
                      soesales + relationship + licenses + interaction+
                      TrustCourts * lintensity_investment_employment + (1|province),
                      data = mydata)

#------------
MLM_bi <- lmer(letcs ~ 
                 lintensity_investment_employment +
                    marketsize + growthrate + lscale + soe_influence +
                    mgovthelp + mtaxrate  +
                    BI_mc + gdpper_1999_2003 + pop +
                    soe + collective + private + foreign +
                    lage + lemp + sales_otherprov + govtsales + 
                    soesales + relationship + licenses + interaction+
                    BI_mc * lintensity_investment_employment + (1|province),
                  data = mydata)

#------------
MLM_enforce <- lmer(letcs ~ 
                      lintensity_investment_employment +
                 marketsize + growthrate + lscale + soe_influence +
                 mgovthelp + mtaxrate  +
                 enforce + gdpper_1999_2003 + pop + 
                 soe + collective + private + foreign +
                 lage + lemp + sales_otherprov + govtsales + 
                 soesales + relationship + licenses + interaction+
                 enforce * lintensity_investment_employment + (1|province),
               data = mydata)

stargazer(MLM_LSC,MLM_court,MLM_bi, MLM_enforce, keep.stat = ("n"), 
          type = "html", out = "Table3.doc",
          covariate.labels =c("Fixed-Asset Intensity (log)","Market Size","Growth Rate",
                              "Scale Economies (log)","SOE Dominance", "Government Help","Tax Burden",
                              "LSC","Trust in Courts","Bureaucratic Integration", 
                              "Enforcement (Factor Scores)",
                              "GDP per Capita (log)","Population",
                              "SOE","Collective","Private","Foreign","Age (log)",
                              "Employees (log)","Out-of-Province Sales","Proportion of Sales to Government",
                              "Proportion of Sales to SOEs","Years of Relationship",
                              "Licenses","Interaction with Government", 
                              "FAI x LSC", "FAI x Trust in Courts",
                              "FAI x Bureaucratic Integration", "FAI x Enforcement (Factor Scores)"),
          no.space = TRUE)


#===============================================================================
# Figure 1: Plot Marginal Effects of FAI
#===============================================================================
#install.packages("interplot")
#install.packages("gridExtra")

library("interplot")
library("gridExtra")

set.seed(11321)
plot_lsc <- 
  interplot(MLM_LSC,"lintensity_investment_employment","LSC") +
  geom_hline(yintercept = 0,color="white",size=2,alpha=0.7) +
  geom_line(size=1) +
  labs(x="Legal Service Capacity (Mean Centered)",y="Marginal Effects") +
  theme(
    title = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    axis.text.x = element_text(size = 10),
    axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0),size=12),
    axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0),size=12))+
  scale_x_continuous(breaks = round(seq(min(mydata["LSC"]), max(mydata["LSC"]), by = 1),1))

#------------
plot_court <- 
  interplot(MLM_court,"lintensity_investment_employment","TrustCourts") +
  geom_hline(yintercept = 0,color="white",size=2,alpha=0.7) +
  geom_line(size=1) +
  labs(x="Trust in Courts (Mean Centered)",y="Marginal Effects")+
  theme(
    title = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    axis.text.x = element_text(size = 10),
    axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0),size=12),
    axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0),size=12))+
  scale_x_continuous(breaks = round(seq(min(mydata["TrustCourts"],na.rm=T), 
                                        max(mydata["TrustCourts"],na.rm=T), by = 0.2),1))

#------------
plot_bi <- 
  interplot(MLM_bi,"lintensity_investment_employment","BI_mc") +
  geom_hline(yintercept = 0,color="white",size=2,alpha=0.7) +
  geom_line(size=1) +
  labs(x="Bureaucratic Integration (Mean Centered)",y="Marginal Effects")+
  theme(
    title = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    axis.text.x = element_text(size = 10),
    axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0),size=12),
    axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0),size=12))+
  scale_x_continuous(breaks = round(seq(min(mydata["BI_mc"],na.rm=T), 
                                        max(mydata["BI_mc"],na.rm=T), by = .5),1))

#------------
plot_enforce <- 
  interplot(MLM_enforce,"lintensity_investment_employment","enforce") +
  geom_hline(yintercept = 0,color="white",size=2,alpha=0.7) +
  geom_line(size=1) +
  labs(x="Enforcement (Factor Scores)",y="Marginal Effects")+
  theme(
    title = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    axis.text.x = element_text(size = 10),
    axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0),size=12),
    axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0),size=12))+
  scale_x_continuous(breaks = round(seq(min(mydata["enforce"],na.rm=T), 
                                        max(mydata["enforce"],na.rm=T), by = 1),1))

plots.combined<-grid.arrange(plot_lsc, plot_court, plot_bi, plot_enforce , 
             ncol = 2, nrow = 2)

ggsave("marginaleffects.jpg", plot = plots.combined,
       device = "jpg",width=8,height=5.5, dpi = 600)

